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SHOCK CAPTURING SCHEMES FOR MULTIDIMENSIONAL FLOW 


ABSTRACT 

We report on progress made in the development of a ‘genuinely multi- 
dimensional’ finite-volume algorithm for problems in inviscid gasdynamics. 
The approach entails (a) the reconstruction of flowfield data using a planar 
wave pattern in which the strengths and orientations of the component waves 
are derived independently of the mesh geometry, and (b) the development 
of a flux formula which provides a numerical approximation to the flux at 
a finite-volume cell face during the passage of waves which are in general 
oblique to the face. We outline several algorithms, and include the most 
recent developments. The results of several numerical test cases are also 
included. 


I. INTRODUCTION 

The standard approach to modeling Euler flowfields in several space di- 
mensions is to approximate the multidimensional solution operator by a se- 
quence of one-dimensional operators applied along the mesh cooordinate di- 
rections. Several methods of this type are highly developed and have been 
applied successfully to a vast body of flow problems governed by the Euler 
and Navier-Stokes equations. The success of these methods is not, however, 
complete. Solution quality can be, in some instances, strongly mesh depen- 
dent. This is especially true in discontinuous regions whenever the dominant 
wave transition is oblique to the grid. This observation defines an obvious 
goal: to develop a grid-independent high-resolution flow solution algorithm. 

The research undertaken during the grant period and described in this 
report is a preliminary effort in the development of such a scheme. The 
approach developed here is based on a local flowfield reconstruction through 
a pattern of grid-oblique waves. This approach reduces, in a natural way, the 
dependence of the solution on the grid geometry, and represents a significant 
step toward the realization of a grid-independent algorithm. 

Research which predates the present work to improve wave resolution 
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in multidimensional flow includes the work of Davis, 1 Roe, 2 Hirsch et a J., 3 
Powell and Van Leer, 4 and Levy et a J. s These methods are outlined briefly 
below. 
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In the approach suggested by Davis, the passage of waves at an interface 
is modeled in a coordinate system which is locally rotated in such a way that 
the flow is nearly one-dimensional in the rotated frame. The input states 
for this rotationally-biased algorithm are determined by interpolation in the 
vicinity of the cell face. Davis assumes that the primary wave is a shock 
with normal aligned with the velocity jump across the interface. An upwind 
algorithm is used to compute the waves in the primary direction; terms 
which arise in the direction perpendicular to the primary waves are centrally 
differenced. Convergence is enhanced (at the expense of wave resolution) 
by smoothing the wave orientation information, and also freezing the wave 
angles at some stage during the calulations. This method produces solutions 
in which oblique shocks are resolved very well. Davis’s procedure does not, 
however, explicitly account for a single wave of the stream type, because such 
a wave would be perpendicular to the specified shock direction. 

Levy et a 1. have used a similar procedure to compute simple two-dimensional 
flows, and have explored a variety of choices for the local angle of rotation. 

In the continuation of this work, the use of upwinding algorithms in both the 
primary and perpendicular wave directions has been explored. 

Roe’s method entails the reconstruction of gradient data on a triangle. 
Because we use a similar reconstruction technique in one of the algorithms 
developed here, we defer the description of this method to Section II. 


A conservative scheme based on the characteristic form of the Euler equa- 
tions has been proposed by Hirsch et a I. In this scheme, the underlying 
physics is modeled by the characteristic compatibility conditions written for 
the choice of wave normals which minimizes the coupling between waves. In 
other words, the governing system of equations is written in the form 


dW , dW , dW 

"aT + A * aT + A " a7 


= R, 


where W are the characteristic variables, such that R is minimized. Because 
the choice of normals determines the specific waves which are computed, 
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and because this choice depends on the gradients of the flow properties, 
this is a genuinely multidimensional approach. Powell and Van Leer have 
implemented a conservative algorithm based on this technique. 

Other multidimensional algorithms, related to one or more of the methods 
described above, have been recently reported by Dadone and Grossman , 6 and 
by Kontinos and McRae . 7 

We begin with an outline of local flow reconstruction techniques investi- 
gated during the couse of this work. The flux formula is then developed in 
Section III, and the application of the most recently developed algorithm to 
several planar flow test cases is summarized in Section IV. 

II. FLOWFIELD RECONSTRUCTION 

In the finite- volume discretization of the flow problem, the physical region 
of interest is divided into a number of cells, and the distribution of flow 
properties within each cell is then specified according to some convenient 
approximation (for instance, that the flow properties are constant, or vary 
linearly in the cell). The wave system which governs the evolution of the 
flowfield is then constructed by matching the wave transitions to the local 
variations in the flow properties. The flowfield then evolves in a manner 
which is consistent with the passage of these waves, through the influence of 
the waves on the flux at the cell faces. 

In this section we describe two methods by which the flowfield is recon- 
structed locally. The first is a Riemann method which is a generalization of 
the classical grid-aligned wave model. 


Two-State Reconstruction 

Given any two neighboring states in a flowfield, there are an infinite 
number of wave patterns by which the jump in state variables might be 
reconstructed. Clearly, additional rules are necessary to make the choice 
of waves unique. In the classical wave model, the waves are chosen to be 
oriented in a specified (grid-dependent) direction, and this choice, together 
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with a further restriction on the number of waves, leaves precisely enough 
free wave parameters (the wave strengths) to make the wave pattern unique 
for a given jump in states. The primary advantage of this technique lies 
in its simplicity. There is, as pointed out in the Introduction, a drawback 
in this reconstruction technique: High resolution of wave structures can be 
expected only if the true wavefront is aligned with the assumed (grid-aligned) 
wave direction. The failure of a grid-aligned wave model to properly describe 
oblique waves in multidimensional data has been recognized for some time 
(see, for example, Ref.2). 

The investigation of two-state reconstruction techniques undertaken here, 
and reported in Refs. 8 and 9 evolved in parallel with the method of Rum- 
sey et al . 10,11 As an aid to the description of the method, we introduce a 
remarkably useful picture of linearized state space, ^ due to Roe. 12 An arbi- 
trary transition (Su,Sv, Sp) between states / and r in planar flow is depicted 
in Figure 1 in a coordinate system in which the coordinate directions are 
chosen to be the pressure and the x— and y— velocity components. The state 
/ is placed at the origin, and the right state r is at the point (Su, Sv, Sp). The 
envelope of possible transitions through acoustic waves with one endstate at 
/ must lie on a cone with apex at /, because all such transitions satisfy the 
relation 

\Sp\ = />a^/(6u) 2 + (6t>) 2 . (1) 

All transitions across shear waves lie in horizontal planes. Note that den- 
sity transitions are not representable on this diagram; thus an entropy wave 
cannot be depicted. 


A superposition of k elementary wave transitions is used to reconstruct 
an arbitrary jump, that is 

SV = tT6V,. 

J = 1 


The problem is, of course, to determine the number, orientations, and strengths 
of the waves to be used. Assuming that the overriding requirement is single 

^ An efficient approach to shock capturing relies on the use of a local linearization 
technique which retains the important features of the nonlinear problem. 9,10 In what 
follows, we assume a local linearization, so that the state variables are replaced whenever 
appropriate by local averages. (For instance, p and a in (1) are average values). 
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Figure 1: Elementary wave transitions in state space. 

Rankine-Hugoniot wave recognition, any candidate wave model must return 
a single wave for this special case. 

The type of wave pattern used in Refs. 8 - 11 to construct the state jump 
was guided by some measure of the closeness of the union of wave paths to 
the overall jump. If the number of waves is restricted to four (as in the grid- 
aligned case), the minimum path length choice leads to two acoustics and an 
entropy wave for r inside the acoustic cone; for r outside the cone, the waves 
are a forward- or a backward-facing acoustic, and a shear. The wave pattern 
specified here is shown in Figure 2. Note that the velocity jump across each 
acoustic and shear is parallel to the overall velocity jump (and the acoustic 
and shear waves are perpendicular to one another). A variety of other rules 
have been devised, and some of these are reported in Refs. 9-11. 


Three-State Reconstruction 

As we have seen above, two-state methods provide a local description 
of the flowfield based on the jump information between neighboring states. 
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Figure 2: Minimum pathlength wave transitions. 

Our experience with these methods has led to the conclusion that it is un- 
reasonable to expect a realistic description of the local wave content of the 
flowfield based on this very restricted information (test results followed by 
further discussion of these methods appear in the remainder of this report). 

A natural extension of the reconstruction technique entails the introduc- 
tion of a third data point into the local description of the flow; this provides 
sufficient information to approximate flow gradients* This technique was 
first proposed by Roe . 2 

Three-point reconstruction fits naturally into the standard unstructured 
mesh discretization, in which the mesh consists of a set of nonoverlapping 
triangles. The data is assumed to be stored at the nodes; the data provides 
estimates of the flowfield gradients to which are matched the variations in 
flow properties described by the chosen wave pattern. There remains, of 

*The reconstruction of gradients makes it impossible to properly fit discontinuities 
into the data. This drawback is, however, of little importance, because the two-state 
reconstruction procedure is part of a solution algorithm which does not preserve perfectly- 
resolved discontinuities. 
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course, the difficulty of narrowing the choice of the wave pattern to a set 
which is unique for given data, and this issue is central to the reconstruction 
procedure. 




In planar flow, there are eight derivatives, from which eight wave param- 
eters can be found. One might therefore suppose that the natural choice in 
planar flow is a set of four waves (each of unknown strength and orientation), 
but a study of the planar-wave solutions of the Euler equations shows that 
the minimum number of waves (for arbitrary data) is in fact five, consisting 
of two acoustics, two shear waves, and an entropy wave. This is the choice 
we make in the algorithm described in Ref. 13. It is possible, of course, to 
specify a wave pattern with more than five waves; Roe and co-workers have 
proposed several six-wave models. 2 ’ 14 Finally, we observe that a reasonable 
choice of waves must also lead to simple algebra, otherwise the efficiency of 
the algorithm will suffer. 


Because only eight wave parameters can be specified independently, addi- 
tional rules are necessary to reduce the number of unknown wave parameters 
to precisely eight. The five-wave model proposed in Ref. 13 consists of a 
forward- and backward-facing acoustic wave pair with common normal, and 
a mutually perpendicular shear wave pair. This choice leads to the wave 
normals 


n 


AC 



and 


the first shear wave angle is 


Vp - p/(7p)Vp 

I Vp - p/( 7 p) Vp|’ 


^shearj — ~ <H*CtcLIl 


{Vy ^ar) COS (u x "I" Vy) 
(Vj- H” Wy) “f" Sill 20ac {\tx + Vy) 


It is possible to derive expressions for the wave strengths directly from the 
data on the triangle, but a more robust algorithm results if the triangle data is 
now reconstructed by an edge-by-edge procedure, in which the wave patterns 
between pairs of states consist of waves with the above-specified orientations, 
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but with strengths found from the edge state jump. The expressions for the 
wave strengths a in terms of the state transition SV along an edge are 

= £ ftp ± pa 6u ■ n“] , 

<*en = Sp- —tip, 

7P 


and 

where 


a,hi = -n ,h * • 8u', = n* hl • 6u', 

8u' — 8u — pa (a«c, — a«c a ) w* 0 


This edgewise procedure leads to an inconsistency, because the average 
states on the edges differ from the triangle-averaged state. This issue remains 
unresolved. 


III. FLUX FORMULA 

In this section, we give an account of the development of formulas to 
predict the numerical flux function. We begin with the flux function used in 
conjuction with the two-state reconstruction described earlier in this report. 

Figure 3 depicts adjacent cells in the mesh; the flux is to be estimated 
at the common face. Suppose the wave pattern consists of mutually perpen- 
dicular waves, and we denote the orthogonal wave-aligned directions by the 
unit vectors i n and i t . The upwind flux formula used in Refs. 8-11 can be 
interpreted as the flux estimate which arises from the assumption that the 
flux components in the wave aligned directions are 

F" = I [F,“ + F„" - £ |*f;i] , 

where the sum is over waves with normal t n , and 

F' = \ [F! + F! - £ |«Fi|] , 
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Figure 3: Figure for the derivation of the flux formula (2) 

where the sum is over waves with normal it- The face- normal flux is then 

F 1 = F*i t • n } + F n i n -n s = \ [f/ + F/ - £ \6F w i w ■ n f \] , (2) 

where n } is the cell face normal. This flux formula has also been applied in 
algorithms in which the two-state reconstruction includes waves of arbitrary 
relative orientation (for which the proper interpretation of the formula is not 
clear). As we shown in Section IV below, the use of this formula leads, in 
general, to highly nonmonotonic solutions. 

A different flux formula, which has a clearer interpretation for waves of 
arbitrary orientation, has been developed for use in conjuction with the three- 
state reconstruction describe in Section III. The details of the development 
we outline below are presented in Ref. 13. The reader is also referred to 
Roe’s original paper. 2 

We assume that the finite- volume cells are built from the triangle-centroid 
mesh dual shown in Figure 4. Then the cell faces are the edges of the polyg- 
onal cells, and we wish to estimate the flux on each face. Because the wave 
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mesh dual 



Figure 4: Finite- volume cells. 

reconstruction is applied to the data on a triangle, we divide each cell face 
into the two pieces that are contained in adjacent triangles. The flux for- 
mula is then used to estimate the flux on partial-faces (which is equivalent 
to writing the flux on the full face as a weighted arithmetic average of the 
flux predicted on the two partial-faces). 

Suppose each finite volume cell in Figure 4 is assembled from a set of 
triangles which are constructed by connecting the vertices of the polygonal 
cell to the center node. Such a triangle is shown in Figure 5. We now imagine 
a single plane wavefront passing through this triangle. Let the state ahead 
of the wave be denoted by the subscript r, and that behind the wave by /. 
The contour integral of the normal component of the flux tensor around the 
cell gives the simple result 

j t [Qa) = -s w 6F”, (3) 

where A is the area of the triangle, s w is the length of the wavefront inside 
the triangle, 8F W is the normal flux jump across the wave, and Q is the 
area-averaged value of the vector of conserved variables (mass, momentum, 
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Figure 5: Wavefront passing through a triangular cell fragment 

and energy per unit area). The result (3) is used below in the development 
of the numerical flux function. 

Figure 6 depicts two triangular fragments of neighboring cells l and r. 
For algorithmic simplicity, we wish to account for the entire effect of the 
waves which arise between the states l and r through the numerical flux at 
the common face. This would obviate the need to explicitly calculate the 
effect of these waves on the remaining faces of the triangular segments. The 
manner in which this can be accomplished is to set the flux at the common 
face, normal to the face, to 

s ‘ 

F = F t + 6F W —, (4) 

s f 

where s l w is an estimate of the length of the wave segment contained in the 
left triangle, and sj is the face length. We note that it is equally valid to 
write 

F = F t - 6F w —, 

s f 

and it is perhaps most reasonable to use the average of (4) and (5). Finally, 
for simplicity, we assume that ratio of the length of the wave segment to the 
cell face length is zero or unity, depending on the direction of motion of the 
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Figure 6: Neighboring triangular cell fragments 


wave. This yields the flux formula 





* = 1 


( 6 ) 


where Ajt is the wave speed of the k — th wave, for a system of five waves. 


It is interesting to note that the only difference in the flux formula used 
in Refs. 8-11, which was derived earlier in this section, and can be written 
in the form 




F, + Fr-± sgn(A t )^ | n 


k=l 



and equation (6) is the dot product term. 


IV. TEST CASES 

The methods described in this report have been tested on a array of 
simple two-dimensional flow problems. We show results for three channel 
flows, which cover the Mach number range from subsonic to high-supersonic. 

Simple forward time integration with local time stepping is used in every 
case we show below. 




Figure 7: Channel geometry and grid. 

The two-state reconstruction algorithm tested here on the supersonic 
channel flow differs slightly from earlier implementations: the method is 
applied here on an unstructured grid; the calculations presented in Refs. 8 
and 11 we done on a structured grid. We also note that the finite- volume 
discretization is, for this algorithm, cell-centered. 


Supersonic Channel Flow 

This is the problem of Levy et a J. 5 We present the geometry of the channel 
and a uniform 1074 node grid in Figure 7. The ramp angle on the lower 
wall is 15°. The Mach number contours shown in Figure 8 for an inflow 
Mach number of 2 show the solution obtained using the minimum pathlength 
two-state reconstruction and the flux formula (2). As this figure shows, 
although the wave resolution is good, the solution is clearly nonmonotonic. 
Furthermore, the maximum density residual diminshes only one order of 
magnitude over 1000 timesteps. 

The same case was computed using the three-state reconstruction algo- 
rithm, together with the flux formula (6). Pressure and Mach number con- 
tours in Figures 9 and 10. These figure show that the variations in the flow 
properties are very nearly monotonic, and the resolution of the shockwaves 
is very high (2-3 triangles). The density residual history, shown in Figure 
11, indicates better than two order-of-magnitude convergence in about 1000 
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Figure 8: Mach number contours. 
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Figure 9: Pressure contours, method of Ref. 13. 
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Figure 10: Mach number contours, method of Ref. 13. 
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Figure 11: Density residual history for the supersonic channel flow problem 
steps. 


Double Ramp 

The monotonicity of the algorithm of Ref. 13 is tested far more severely 
in this flowfield, which is a Mach 4 flow over a 20-35° double ramp. The 
uniform 1052 node grid we used is shown in Figure 12. Pressure and Mach 
number contours are plotted in Figures 13 and 14. The Mach number 
contours in Figure 14 clearly indicate the presence of the slip surface. These 
contour plots also show nearly-monotone variations in the flow properties. 
The density residual history in Figure 15 shows a two order-of-magnitude 
reduction in 500 timesteps. 


Subsonic Channel Flow 

The last test case we show here is a fully subsonic flow in the channel 
shown in Figure 16. The bump on the lower wall is a 10% circular arc, and 


Km QUAU'iY 
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Figure 13: Pressure contours 
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Figure 15 : Density residual history for the double-ramp problem 
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Figure 16: Channel geometry and grid 
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Figure 17: Pressure contours. 

there are 29 nodes on the bump. The stretched grid contains a total of 1108 
nodes. Pressure and Mach number contours for an inflow Mach number of 
0.6 are shown in Figures 17 and 18. The more serious convergence problem 
alluded to earlier is apparent from the density residual history plot in Figure 
19. 


Discussion 

It is immediately clear from these examples that the genuinely multidi- 
mensional methodology leads to very high wave resolution. Unfortunately, 
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Figure 18: Mach number contours. 




Figure 19: Density residual history for the subsonic channel flow problem 
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this comes at the cost of monotonicity and convergence. This is especially 
true of results using the two-state reconstruction method together with the 
flux formula (2). (A detailed study of the monotonicity behaviour of this 
family of methods is presented by Rumsey et al . n ) 

There is a vast improvement in the quality of the solution when a three- 
state reconstruction is used with the new flux formula (6), but the solution 
still fails to be strictly monotonic, and convergence is not enhanced suffi- 
ciently. The lack of monotonicity and convergence is the most apparent in 
subsonic flow, perhaps because of the omnidirectional nature of propaga- 
tion of acoustic signals. The wave model is designed to identify dominant 
waves in the data, and an added level of sophistication may be necessary to 
reconstruct data in which there are no strongly-preferred directions. 

V. CONCLUDING REMARKS 

Two major issues that have been raised repeatedly in the development 
of genuinely multidimensional schemes have yet to be resolved satisfactorily. 
The first of these is the nonmonotone behavior of solutions, and the second 
is the failure of the solution to converge to steady state. The developments 
described in this report address both of these issues, but it is clear that much 
remains to be done before methods of this type become viable replacements 
for classical grid-aligned schemes. On the other hand, genuinely multidimen- 
sional algorithms (both the finite- volume schemes discussed here, and the 
multidimensional fluctuation distribution algorithms proposed by Roe and 
co-workers) do show the potential to dramatically improve wave resolution, 
and at very little cost. 
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